In [1]:
import os
%matplotlib inline
from netCDF4 import Dataset
from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
In [2]:
#gunung = 'gn_anak_krakatau' #(tgl 3 januari)(10juni)
#gunung1 = 'Gunung Krakatau'
#lat1, lat2 = 1150, 1250
#lon1, lon2 = 1200, 1350
#x, y = 98.392, 3.17
#gunung = 'gn_agung' #(tgl 12-13,17-18,24 mei)(10juni)
#gunung1 = 'Gunung Agung'
#lat1, lat2 = 1050, 1150
#lon1, lon2 = 1700, 1900
#lon1, lon2 = 1700, 1850
#x, y = 115.503622, -8.340686
gunung = 'gn_sinabung' #(7 mei),(9 juni)
gunung1 = 'Gunung Sinabung'
lat1, lat2 = 1600, 1700
lon1, lon2 = 850, 1000
x, y = 98.392, 3.17
#gunung = 'gn_dukono' #(9-10,16-17,23-25 mei)
#gunung1 = 'Gunung Dukono'
#lat1, lat2 = 1500, 1650
#lon1, lon2 = 2300, 2500
#x, y = 127.88, 1.68
#directory = 'C:/Users/Szczynk/Desktop/PKL/data NC 9 juni studi kasus gn.sinabung/'
#img_dir = 'C:/Users/sugeng/Downloads/Bagus Nuryasin/RST ash TIR dan MIRTIR Gunung Dukono/'
#mean_dir = 'C:/Users/sugeng/Downloads/Bagus Nuryasin/'
img_dir = 'C:/Users/Szczynk/Desktop/Bagus Nuryasin (2)/RST ash TIR dan MIRTIR '+gunung1+'/'
mean_dir = 'H:/'
jam = ['0000','0010','0020','0030','0040','0050','0100','0110','0120','0130','0140','0150','0200','0210','0220','0230','0240','0250','0300','0310','0320','0330','0340','0350','0400','0410','0420','0430','0440','0450','0500','0510','0520','0530','0540','0550','0600','0610','0620','0630','0640','0650','0700','0710','0720','0730','0740','0750','0800','0810','0820','0830','0840','0850','0900','0910','0920','0930','0940','0950','1000','1010','1020','1030','1040','1050','1100','1110','1120','1130','1140','1150','1200','1210','1220','1230','1240','1250','1300','1310','1320','1330','1340','1350','1400','1410','1420','1430','1440','1450','1500','1510','1520','1530','1540','1550','1600','1610','1620','1630','1640','1650','1700','1710','1720','1730','1740','1750','1800','1810','1820','1830','1840','1850','1900','1910','1920','1930','1940','1950','2000','2010','2020','2030','2040','2050','2100','2110','2120','2130','2140','2150','2200','2210','2220','2230','2240','2250','2300','2310','2320','2330','2340','2350']
bln = '06'
tgl = ['01','02','03','04','05','06','07','08','09','10','11','12','13','14','15','16','17','18','19','20','21']
thn = '2019'
print(jam.index('1500'))
print(len(tgl[8:9]))
print(tgl[8],bln,thn)
In [3]:
#img_r = np.load(mean_dir+'data_nc_'+gunung+'_B13-B15_1-17_'+bln+'_2019.npy')
#img_g = np.load(mean_dir+'data_nc_'+gunung+'_B07-B13_1-17_'+bln+'_2019.npy')
img_r = np.load(mean_dir+'data_nc_'+gunung+'_B13-B15_1-17_05_2019.npy')
img_g = np.load(mean_dir+'data_nc_'+gunung+'_B07-B13_1-17_05_2019.npy')
mean_R = np.mean(img_r, axis=0)
std_R = np.std(img_r, axis=0)
mean_G = np.mean(img_g, axis=0)
std_G = np.std(img_g, axis=0)
#low, mid, high, extreme = -0.5, -1.5, -2.5, -3.5
low, mid, high, extreme = 0., -1., -2., -3.
tgl1 = 8
for l in range(len(tgl[8:9])):
for n in range(len(jam)):
try:
directory = 'H:/dataNC/'+bln+'/'+tgl[tgl1+l]+'/'
dsetB07 = Dataset(directory+'H08_B07_Indonesia_'+thn+bln+tgl[tgl1+l]+jam[n]+'.nc', mode = 'r')
dsetB13 = Dataset(directory+'H08_B13_Indonesia_'+thn+bln+tgl[tgl1+l]+jam[n]+'.nc', mode = 'r')
dsetB15 = Dataset(directory+'H08_B15_Indonesia_'+thn+bln+tgl[tgl1+l]+jam[n]+'.nc', mode = 'r')
lat = dsetB13.variables['latitude']
lon = dsetB13.variables['longitude']
B07 = dsetB07.variables['I4'][...,lat1:lat2,lon1:lon2]
B13 = dsetB13.variables['IR'][...,lat1:lat2,lon1:lon2]
B15 = dsetB15.variables['I2'][...,lat1:lat2,lon1:lon2]
B07 = np.squeeze(B07)
B13 = np.squeeze(B13)
B15 = np.squeeze(B15)
R = B13.data - B15.data
G = B07.data - B13.data
r = R
TIR = (r-mean_R)/std_R
#TIR = ndimage.median_filter(TIR, size=3)
g = G
MIRTIR = (g-mean_G)/std_G
res = np.zeros((len(B13), len(B13[0])))
for i in range(len(B13)):
for j in range(len(B13[0])):
if TIR[i,j] < high and MIRTIR[i,j] > 0.:
res[i,j] = 3.
elif TIR[i,j] < mid and MIRTIR[i,j] > 0.:
res[i,j] = 2.
elif TIR[i,j] < low and MIRTIR[i,j] > 0.:
res[i,j] = 1.
else:
res[i,j] = 0.
m = Basemap(projection='merc', resolution='i',llcrnrlon=lon[lon1], llcrnrlat=lat[lat1], urcrnrlon=lon[lon2] ,urcrnrlat=lat[lat2])
x1 = np.linspace(lon[lon1],lon[lon2],len(B13[0]))
y1 = np.linspace(lat[lat1],lat[lat2],len(B13))
xx, yy = np.meshgrid(x1, y1)
x2, y2 = m(xx,yy)
#lev = [0,1,2,3,4]
lev = [0,0.99,1.99,2.99,3.99]
colva = [(1,1,1),(1,1,0),(1,0.7,0),(1,0,0)]
plt.figure(figsize=(11,11))
m.drawcoastlines()
m.drawcountries()
debu = m.contourf(x2, y2, res, levels=lev, colors=colva)
cbar = m.colorbar(debu)
cbar.set_ticklabels(['Tidak','Rendah', 'Menengah', 'Tinggi'])
locgn = m.scatter(x, y, marker='^', s=10, color='m', latlon=True)
plt.legend([locgn],[gunung1])
plt.title(gunung1+' '+tgl[tgl1+l]+'-'+bln+'-'+thn+'-'+jam[n]+' UTC')
plt.savefig(img_dir+'RST ash '+gunung+' '+tgl[tgl1+l]+'-'+bln+'-'+thn+'-'+jam[n]+' UTC.png',bbox_inches='tight')
except:
pass
print(tgl[tgl1+l])
In [4]:
m.drawcoastlines()
m.drawcountries()
locgn = m.scatter(x, y, marker='^', s=10, color='m', latlon=True)
plt.legend([locgn],[gunung1])
plt.title('Daerah '+gunung1)
plt.show()
In [5]:
#plt.figure(figsize=(13,14))
plt.subplot(221)
m.drawcoastlines()
m.drawcountries()
locgn = m.scatter(x, y, marker='^', s=10, color='m', latlon=True)
plt.legend([locgn],[gunung1])
m.imshow(mean_R[:,:], cmap='coolwarm')
m.colorbar()
plt.title('Mean R')
plt.subplot(222)
m.drawcoastlines()
m.drawcountries()
locgn = m.scatter(x, y, marker='^', s=10, color='m', latlon=True)
plt.legend([locgn],[gunung1])
m.imshow(mean_G[:,:], cmap='coolwarm')
m.colorbar()
plt.title('Mean G')
plt.subplot(223)
m.drawcoastlines()
m.drawcountries()
locgn = m.scatter(x, y, marker='^', s=10, color='m', latlon=True)
plt.legend([locgn],[gunung1])
m.imshow(std_R[:,:], cmap='coolwarm')
m.colorbar()
plt.title('Std R')
plt.subplot(224)
m.drawcoastlines()
m.drawcountries()
locgn = m.scatter(x, y, marker='^', s=10, color='m', latlon=True)
plt.legend([locgn],[gunung1])
m.imshow(std_G[:,:], cmap='coolwarm')
m.colorbar()
plt.title('Std G')
plt.show()
In [ ]: